""" Our interpolation functions

    :Authors:
      - 20081016 Nicola Creati
      - 20081016 Roberto Vidmar

    :Revision:  $Revision: 52 $

    :Copyright: 2011-2012
                Nicola Creati <ncreati@inogs.it>
                Roberto Vidmar <rvidmar@inogs.it>

    :License: MIT/X11 License (see :download:`license.txt
                               <../../license.txt>`)
"""

import numpy as np

def interp2d(z, ext, xy):
  """ Return a new 2d matrix of linear inetrpolated values at xy[i, j].

      z is a 2d array of values at evenly spaced values in
      ext = ((xmin, xmax), (ymin, ymax))

      :param z: 2d array of evenly spaced values
      :type z: numpy array
      :param ext: ((xmin, xmax), (ymin, ymax))
      :type ext: tuple
      :param xy: 2d matrix
      :type xy: numpy array
      :returns: new 2d matrix of linear inetrpolated values at xy[i, j]
      :rtype: numpy array
      :raises:

      .. warning:: All xy points must fit into ext rectangle.
                   **NO CHECK IS MADE!**
  """
  (xmin, xmax), (ymin, ymax) = ext
  if xmin > xmax or ymin > ymax:
    raise ValueError ("Wrong limits: %s > %s or %s > %s" % (
      xmin, xmax, ymin, ymax))
  nrows, ncols = z.shape

  # Compute grid spacing in both directions
  dx = (xmax - xmin) / (ncols - 1.)
  dy = (ymax - ymin) / (nrows - 1.)

  c, t = divmod((xy[:, 0] - xmin), dx)
  r, u = divmod((xy[:, 1] - ymin), dy)

  t = t / 10.
  u = u / 10.

  r = r.astype(int)
  c = c.astype(int)

  newz      = ((1 - t) * (1 - u) * z[r,     c    ] +
                  t    * (1 - u) * z[r,     c + 1] +
                  t    *    u    * z[r + 1, c    ] +
               (1 - t) *    u    * z[r + 1, c + 1])
  return  newz

if __name__ == '__main__':

  z = np.array([[0.0, 1.0, 2.0, 3.0, 4.0],
                [0.0, 1.0, 2.0, 3.0, 4.0],
                [0.0, 1.0, 2.0, 3.0, 4.0],
                [0.0, 1.0, 2.0, 3.0, 4.0],
                [0.0, 1.0, 2.0, 3.0, 4.0]])

  zext = ((10, 20), (30, 40))
  print zext
  xy = np.array([[15.0, 35.0], [12., 31.], [18.5, 39.5], [16.45, 36.1]])
  print "xy=", xy
  print "z=", z

  newz = interp2d(z, zext, xy)
  print "newz=", newz
